---
title: "Humidity interpolation"
author: "Francesco Bailo"
date: "08/07/2020"
output: html_document
---

```{r}
library(raster)
library(sp)
library(sf)
library(meteoland)
library(dplyr)
library(biogeo)
library(ggplot2)
library(gridExtra)
```


# SpatialGridDataFrame

```{r}
ita.sp <- 
  readRDS("gadm36_ITA_0_sp.rds")

dem.r <- 
  raster("30n000e_20101117_gmted_med300.tif")

dem_ita.r <- 
  crop(dem.r, ita.sp)

dem_ita.sp <- 
  as(dem_ita.r, 'SpatialGridDataFrame')

summary(dem_ita.sp)
```

# Station Data

```{r}
load("meteonetwork/station_details.RData")
load("meteonetwork/station_rank.RData")

station_rank$`Anno Solare` <- as.numeric(gsub("\\(|\\)|%", "", station_rank$`Anno Solare`)) 
station_rank$`Ultimo mese` <- as.numeric(gsub("\\(|\\)|%", "", station_rank$`Ultimo mese`)) 

top_quality_stations <- 
  station_rank$Luogo[station_rank$`Anno Solare` >= 95 &
                       station_rank$`Ultimo mese` >= 95]

station_details$lon <- as.numeric(gsub(" E", "", station_details$Longitudine))
station_details$lat <- as.numeric(gsub(" N", "", station_details$Latitudine)) 
station_details$elevation <- as.numeric(gsub(" mslm", "", station_details$Altitudine))

ggplot() + 
  geom_sf(data = st_as_sf(ita.sp), fill = NA) +
  geom_point(data = station_details %>% 
               filter(code %in% top_quality_stations),
             aes(x=lon,y=lat))
```

```{r}
load("meteonetwork/weather_data.df.RData")

weather_data.df <- 
  weather_data.df %>%
  dplyr::filter(code %in% top_quality_stations) %>%
  dplyr::mutate(date = as.Date(date, format = "%d-%m-%Y"), 
                tmin = as.numeric(gsub("°C", "", `Temp. Min`)),
                tmax = as.numeric(gsub("°C", "", `Temp. Max`)),
                prec = as.numeric(gsub(" mm/h", "", Pioggia)),
                hmin = as.numeric(gsub("%", "", `Umid. Min`)),
                hmax = as.numeric(gsub("%", "", `Umid. Max`)),
                wmax = as.numeric(gsub(" km/h", "", `Raffica Max`)))

# Cleaning
weather_data.df$tmin[weather_data.df$tmin == 0] <- NA
weather_data.df$tmax[weather_data.df$tmax == 0] <- NA
weather_data.df$hmin[weather_data.df$hmin == 0] <- NA
weather_data.df$hmax[weather_data.df$hmax == 0] <- NA
weather_data.df$wmax[weather_data.df$wmax == 0] <- NA

weather_data.df$prec[weather_data.df$hmin == 0 | 
                       weather_data.df$hmax == 0 | 
                       weather_data.df$tmin == 0 |
                       weather_data.df$tmax == 0] <- NA

grid.arrange(
  ggplot(weather_data.df) + geom_density(aes(tmin)),
  ggplot(weather_data.df) + geom_density(aes(tmax)),
  ggplot(weather_data.df) + geom_density(aes(prec)),
  ggplot(weather_data.df) + geom_density(aes(hmin)),
  ggplot(weather_data.df) + geom_density(aes(hmax)),
  ggplot(weather_data.df) + geom_density(aes(wmax)),
  ncol = 2)

weather_data.df <- 
  weather_data.df %>%
  dplyr::group_by(date, code) %>%
  dplyr::mutate(rhum = mean(c(hmin, hmax)))

```

# Model

```{r}
model_stations <- unique(weather_data.df$code)

station_details <- 
  station_details %>% 
  dplyr::filter(code %in% model_stations)

station_details.sp <- 
  SpatialPoints(station_details[,c("lon", "lat")],
                proj4string = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))

tmin.m <- 
  weather_data.df %>%
  dplyr::filter(code %in% station_details$code) %>%
  dplyr::select(code, tmin, date) %>%
  tidyr::pivot_wider(id_cols = code, names_from = date, values_from = tmin) %>%
  dplyr::arrange(match(code, station_details$code)) %>%
  as.data.frame()
row.names(tmin.m) <- tmin.m$code
tmin.m$code <- NULL
tmin.m <- as.matrix(tmin.m)

tmax.m <- 
  weather_data.df %>%
  dplyr::filter(code %in% station_details$code) %>%
  dplyr::select(code, tmax, date) %>%
  tidyr::pivot_wider(id_cols = code, names_from = date, values_from = tmax) %>%
  dplyr::arrange(match(code, station_details$code)) %>%
  as.data.frame()
row.names(tmax.m) <- tmax.m$code
tmax.m$code <- NULL
tmax.m <- as.matrix(tmax.m)

prec.m <- 
  weather_data.df %>%
  dplyr::filter(code %in% station_details$code) %>%
  dplyr::select(code, prec, date) %>%
  tidyr::pivot_wider(id_cols = code, names_from = date, values_from = prec) %>%
  dplyr::arrange(match(code, station_details$code)) %>%
  as.data.frame()
row.names(prec.m) <- prec.m$code
prec.m$code <- NULL
prec.m <- as.matrix(prec.m)

rhum.m <- 
  weather_data.df %>%
  dplyr::filter(code %in% station_details$code) %>%
  dplyr::select(code, rhum, date) %>%
  tidyr::pivot_wider(id_cols = code, names_from = date, values_from = rhum) %>%
  dplyr::arrange(match(code, station_details$code)) %>%
  as.data.frame()
row.names(rhum.m) <- rhum.m$code
rhum.m$code <- NULL
rhum.m <- as.matrix(rhum.m)

wmax.m <- 
  weather_data.df %>%
  dplyr::filter(code %in% station_details$code) %>%
  dplyr::select(code, wmax, date) %>%
  tidyr::pivot_wider(id_cols = code, names_from = date, values_from = wmax) %>%
  dplyr::arrange(match(code, station_details$code)) %>%
  as.data.frame()
row.names(wmax.m) <- wmax.m$code
wmax.m$code <- NULL
wmax.m <- as.matrix(wmax.m)

```

```{r}
interpolator <- 
  MeteorologyInterpolationData(station_details.sp, 
                               elevation = station_details$elevation,
                               MinTemperature = tmin.m,
                               MaxTemperature = tmax.m,
                               Precipitation = prec.m,
                               RelativeHumidity = rhum.m,
                               WindSpeed = wmax.m)
```

```{r}
spatial_coverage <- 
  interpolation.coverage(interpolator, type = 'spatial')
head(spatial_coverage)
```

## Calibration

```{r}
tmin_cal <- 
  interpolation.calibration(interpolator, variable = "Tmin",
                            N_seq = 20,
                            alpha_seq = seq(5, 10, by = 1),
                            verbose = TRUE)

interpolator@params$N_MinTemperature = tmin_cal$N
interpolator@params$alpha_MinTemperature = tmin_cal$alpha

tmax_cal <- 
  interpolation.calibration(interpolator, variable = "Tmax",
                            N_seq = 20,
                            alpha_seq = seq(5, 10, by = 1),
                            verbose = TRUE)

interpolator@params$N_MaxTemperature = tmax_cal$N
interpolator@params$alpha_MaxTemperature = tmax_cal$alpha

save.image()
```

# Interpolation on grid

## Testing

```{r}
library(sf)

prov.sf <- 
  read_sf("../shapefile/Italy/Limiti01012017_g/ProvCM01012017_g/ProvCM01012017_g_WGS84.shp")

dem_no.r <- 
  crop(dem.r,
       as_Spatial(prov.sf %>% 
                    filter(COD_PROV == 3) %>%
                    st_transform(4326)))

dem_no.sp <- 
  as(dem_no.r, 'SpatialGridDataFrame')

sgt <-
  SpatialGridTopography(as(dem_no.sp, "SpatialGrid"), 
                        elevation = dem_no.sp$X30n000e_20101117_gmted_med300,
                        proj4string = dem_no.sp@proj4string)

spplot(sgt, "elevation")

dates <- as.Date(c("2020-02-03", "2020-06-03"))
sum(!(dates %in% interpolator@dates))

ml <- interpolationgrid(interpolator, sgt, dates)

spplot(ml, 2, "MinTemperature")
spplot(ml, 2, "MeanRelativeHumidity")
```


```{r, eval = F, echo = T}

read_chunk('interpolate.R')

```

